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Abstract 

We apply the coupled dynamics of time-dependent density functional theory and Maxwell equa- 
tions to the interaction of intense laser pulses with crystalline silicon. As a function of electromag- 
netic field intensity we see several regions in the response. At the lowest intensities, the pulse is 
reflected and transmitted in accord with the dielectric response, and the characteristics of the en- 
ergy deposition is consistent with two-photon absorption. The absorption process begins to deviate 
from that at laser intensities ~ 10 13 W/cm 2 , where the energy deposited is of the order of 1 eV per 
atom. Changes in the reflectivity are seen as a function of intensity. When it passes a threshold of 
about 3 x 10 12 W/cm 2 , there is a small decrease. At higher intensities, above 2 x 10 13 W/cm 2 , the 
reflectivity increases strongly. This behavior can be understood qualitatively in a model treating 
the excited electron-hole pairs as a plasma. 
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I. INTRODUCTION 



The Maxwell equations describe propagation of electromagnetic fields in bulk matter 
taking into account the material properties by the constitutive relations. For ordinary light 
pulses, the response of the medium is linear in the electromagnetic field and is characterized 
by the linear susceptibilities. In recent experiments with intense and ultrashort laser pulses, 
however, one often encounter conditions which require theoretical treatments beyond the 
linear response. If the perturbative expansion is no longer useful, one need to go back to 
the time-dependent Schrodinger equation for electrons and to solve it in time domain. 

In the last two decades, computational approaches to solve the time-dependent 
Schrodinger equation under intense electric fields have been developed for atoms and small 
molecules [l-3|. For electron dynamics in bulk matter as well as in molecules, one often 
needs to go to the less demanding approach based on time-dependent density-functional 
theory (TDDFT) j^-ll]. We consider the TDDFT is the only ab-initio quantum method 
applicable to high fields in condensed media. 

In this paper, we develop a formalism and computational method to describe propagation 
of intense electromagnetic field in the condensed medium incorporating feedback of electron 
dynamics to the electromagnetic field. This requires a consistent treatment of electrons 
and the electromagnetic field in coupled equations of motion. Such attempts have been 
undertaken by several groups, for isolated molecules [l2|, nano-particles 
16]. 
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151 ] . and gases 



Experimentally, electron-hole plasmas are generated by irradiating so 
pulses, and the threshold for dielectric breakdown has been measured 



ids with strong laser 



17-261. To describe 



the phenomena, mode 
been developed 
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approaches such as a rate equation for electronic excitations have 



32] . For this problem, we have developed a first-principles approach 



331 ]. using 



9), |33| . We calculated dielectric breakdown in crystalline diamond |9| and quartz 
TDDFT and treating the electric field as a longitudinal field. The calculated dielectric 
breakdown threshold was much higher than observed. In the present work, we improve the 
theory by incorporating both magnetic and electric fields in the equations, permitting a 
proper description of transverse electromagnetic wave propagation. 

Our formal development gives a way to separate out the two spatial scales that must be 
treated simultaneously. The electron dynamics is calculated on the atomic scale, resolving 
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position dependences of some tenths of an atomic unit. The electronic field is decomposed 
into two parts, one on the atomic scale and the other on the scale of the electromagnetic 
wave length. The atomic scale field is very similar from one unit cell to neighboring cells of 
the crystal. The other part of the field gives the large scale variation needed to describing 
the electromagnetic self-coupling and wave propagation. We introduce two grid systems 
with different resolution for this problem. 

The construction of the paper is as follows. In Sec. 2, we present our formalism of 
multi-scale description for coupled dynamics of electrons and electromagnetic fields. In 
Sec. 3, numerical methods are explained. In Sec. 4, calculated results are presented. We 
provide an interpretation for the electron dynamics at the surface in terms of the electron- 
hole plasma. We also compare dynamics of the present multi-scale calculation with the 
microscopic dynamics in longitudinal and transverse geometries. Finally a summary is 
presented in Sec. 5. 



II. FORMALISM 



A. Macroscopic equations for electromagnetic fields 

We will consider a coupled dynamics of electrons and electromagnetic field in bulk crys- 
talline solid allowing for strong electromagnetic fields. We immediately recognize there are 
two different spatial scales in the problem. The spatial scale of electromagnetic field is set by 
laser wavelength, of the order 1 /iin. The spatial scale of electron dynamics is much smaller, 
of the order of 1CT 1 nm. We are thus led to a multi-scale description for the problem, em- 
ploying two spatial grids of different grid sizes. We will use the notation R and f for the 
macroscopic and microscopic coordinates, respectively. 

The essence of our method is to use the freedom to choose the electromagnetic gauge to 
separate the two scales. In the expression for the electric field, 

±f (i) 

c at 

the gauge field A contains all the macroscopic electromagnetic physics. The microscopic 

— * 

physics, to be calculated on a unit cell of the lattice, uses both A and the scalar potential <p. 

To derive the theory formally, we start by taking a specific gauge condition, the scalar 
potential <fi is set equal to zero. In this gauge, we have the following equations for the vector 



potential A(f, t), 

1 d - - 

- - ttt V A = Aire (n ion - n) , (2) 
c at 

?flF- v ' A+v ( v - i ) = -— J - (3) 

Here we introduced the ionic density given by 



n : , 



l (r) = J2ZJ(f-R a ), (4) 



where i? a and Z a are the coordinate and charge number of a-th ion, respectively We ignore 
the motion of ions throughout this paper, n and j are the number density and current of 
electrons, respectively, and satisfy the equation of continuity, 

+ Vj = 0. (5) 

We note that longitudinal part of the vector potential is described redundantly by two 
equations (J2J) and (J3J). 

We then proceed to define macroscopic quantities from the microscopic ones. As will be 
discussed later, the microscopic density and current are obtained from the time-dependent 
Kohn-Sham orbitals. The macroscopic version of these quantities, iV R (t) and Jr(£), may 
be defined in principle by applying some smoothing function to the microscopic quantities. 
In practice, we define this by averaging over the unit cell of the lattice. The macroscopic 
density and current also satisfies the equation of continuity. We will not need it for the 
geometry considered below, but it would be needed for other geometries. 

We obtain macroscopic vector potential from A(f, t) by a smoothing procedure which we 
denote as An(t). It satisfies the equations 

--|-V R I R (t) = -47reiV R (t), (6) 

c at 

^Mt) - V R ^ R (t) + V R (V R • A R (t)) = —J*®, (7) 

This is our basic equation to describe a propagation of macroscopic electromagnetic field. 

For a microscopic physics, we treat the macroscopic field as uniform and otherwise we 
retain only a longitudinal part of the vector potential. Physically, the approximation is to 
neglect the transverse current and variation of the magnetic field within the unit cells. The 
neglect of transverse current amounts to ignore the orbital magnetization. The neglect of 



magnetic field effects on electrons may be justified when the velocity of electrons accelerated 
by the laser pulse is much smaller than the velocity of light. As will be explained later, we 
will employ a periodic scalar potential instead of the longitudinal vector potential for the 
microscopic description in the unit cells. 

In general, the presence of boundaries requires special attention. If we may assume that 
a surface charge is localized in a sufficiently thin layer at a surface, we may treat it as a 
discontinuity of the macroscopic vector potential at the surface. Let us consider a small 
volume around a point R at a surface and apply the Gauss theorem to Eq. (jSJ). We obtain 



where n is a unit vector normal to the surface at R and £r(£) is the surface charge at R. 
This equation describes the boundary condition for the macroscopic vector potential across 
the surface. In the geometry we consider here, however, the fields are all parallel to the 
surface so that the macroscopic vector potential is continuous at the surface. 

B. Microscopic equations for electrons 

For the microscopic electron dynamics, we assume a periodic band structure and apply 
the equations of motion of the time-dependent density functional theory. We will make 
several assumptions here. 

First we assume that electron dynamics at different macroscopic positions may be de- 
scribed independently. Namely, we define Kohn-Sham orbitals at every macroscopic grid 
point and ignore any direct interactions between electrons belonging to different macro- 
scopic grid points. We only take into account the interaction between electrons of different 
macroscopic grid points through the macroscopic vector potential A R (£). 

Second, we assume that iV R (i) is independent of time. This condition is satisfied in the 
one-dimensional propagation of linearly polarized light at normal incidence on an interface, 
since the macroscopic current does not include any longitudinal component as discussed 
below. The orbitals evolve under the time-dependent Kohn-Sham equations, but the number 
of electrons in each cell remains the same, and in fact the orbital occupation numbers remain 
zero or one in the time-evolved basis. 

Third, within each cell of the microscopic scale, we ignore any effects of magnetic fields 




(8) 
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on the electrons. The macroscopic vector potential A-&(t) will be treated as a uniform field 
in the microscopic scale. This permits us to treat electron dynamics induced by a uniform 
electric field. We also ignore the transverse component of the microscopic vector potential, 
retaining only the longitudinal part as mentioned before. 

Since all that matters are the physical fields, we are permitted to make a different choice 
of gauge for the microscopic fields. In Eqs. (j2J) and (jHJ), we had chosen the gauge condition 
that removes the scalar potential. However, to take advantage of the periodicity of the 
lattice, we make a gauge transformation at each macroscopic grid point, expressing periodic 
electromagnetic field with a scalar potential instead of the longitudinal part of the mi- 
croscopic vector potential. We denote the scalar potential at macroscopic grid point R as 
<t f R.('f f ,i), to indicate that R will just be a parameter in the equation of motion for <p. 

We denote the Kohn-Sham orbitals at a macroscopic coordinate R as ipi^{f,t). Under 
above conditions and assumptions, the time-dependent Kohn-Sham (TDKS) equation may 
be written as 



ih—iP m (r,t) = j— \ -ihV?+-A K (t)) - eMr,t) + -j^j ^, R (r,t). (9) 

In solving Eq. (Q, the macroscopic coordinate R is treated as a parameter. The Kohn- 
Sham Hamiltonian thus defined is periodic in space and one may introduce Bloch functions 



at each time step, app 



each microscopic cell 



ying periodic boundary conditions on the electron orbitals within 
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34] 



The electron density and current are both periodic in space and are given by 

n R (f,t) = J2\^n(r,t)\ 2 , (10) 

i 

Mr, t) = ^-Y. {^,n(r, t) (-iKV? + ^R (*)) ^,r(^ *) 

-ipi,n{r,t) (W r -+ ^rW) ^r(^*)} > (11) 

where the sum i is over occupied orbitals. The scalar potential 0R,(r, t) satisfies the Poisson 
equation, 

V^> R (f, t) = -4tt (erw )R (r) - en n (r, t)) , (12) 
where rii on ^ is the ionic density at macroscopic grid point R. 
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Since the density and the current are periodic in space, the average over the unit cell is 
meaningful also on the macroscopic scale. The main macroscopic quantity we need from the 
electronic dynamics is the current, defined as 

Jn(t) = ^f n drj R (f,t), (13) 

where Q is the volume of the unit cell. 



C. Conserved energy 



To obtain an expression for the conserved energy in the present multi-scale description, 
we first note the above equations of motion may be derived from the following Lagrangian. 



- df{(en ion - en R ) R - £ IC [n R ]} 



From this Lagrangian, one may derive the following Hamiltonian, 



H 



= fdR Y I df— 



-ihVf + -A R ) 



(14) 



(15) 



The energy calculated from this Hamiltonian is conserved by the equations of motion. 



D. One-dimensional propagation 

In this paper, we will consider a propagation of linearly polarized laser pulse incident 
normally on a bulk crystalline Si on the [110] surface with the laser electric field in [100] 
direction. There are two spatial regions on the macroscopic scale, vacuum and crystalline 
solid. We take a macroscopic coordinate system such that the surface of the crystalline solid 
is xy-plane with z = 0. In this geometry, macroscopic quantities are uniform in both x- 
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and ^-directions. Therefore, macroscopic quantities are specified by the ^-component of R 
which we denote Z. 

The macroscopic vector potential has the following form, 



A R (t) = xA z (t), 



(16) 



In the vacuum region (Z < 0), Az(t) is composed of incident and reflected waves. Inside 
the solid, we assume a locally dipole approximation at each macroscopic coordinate, as 
mentioned before. Then the macroscopic electron current is parallel to the vector potential. 



Mt) = xJ z (t). 



(17) 



We note that this form of electric current is transverse. This justifies our assumption that 
the macroscopic electron density N R (t) is independent of time. 

For later convenience, we summarize equations of motion in the one- dimensional geome- 
try. The vector potential A z (t) satisfies the following equation, 



c 2 dt 2 

The TDKS equation is given by 



1 d 2 . , . d 2 . , . 47re T , . 

Mt) - i&zMt) = —hit)- 



(18) 



d f 1 / - e \ 2 8E 1 

th-^ z (r,t) = I — {-ihV?+ -xA z (t)J - eMr,t) + ^,z(r,t), (19) 



with the density and current, 



n z (r,0=£KMr,0| 2 , 

i 

jz(r, t) = ^- {^*z + -^Azj ^ i)Z - ^,z (ihV? + -^Azj ^*z} , 



■hit) = — J^drxj z (r,t). 



(20) 

(21) 
(22) 



The energy per unit area is a conserved quantity and is given by 

Ea = s / rfz [E /„ *"i | + ^*») *uf 

J^df |-(en ioniZ - en z )0 z + # xc [nz] 



+ 



8vr c V 9t 



<9Z 



(23) 
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E. Linear response 



It is essential that the theory properly describes the propagation of electromagnetic waves 
in the weak field limit, if it is to be useful more generally. In weak fields the electron dynamics 
can be calculated perturbatively to arrive at the usual linear response. The present formalism 
gives the same linear response as other approaches, so the weak field limit will be correct 
if the dielectric function is given correctly. In our previous work, we calculated the linear 
response by separating A into a part that arose from external sources and a part arises from 
the medium [34j. It is not necessary to make this separation in the present formalism. To 
derive the dielectric function, we note Eqs. (IT9"|) . (|2"Tj) . and (I2"2"j) describe relation between the 
macroscopic vector potential Az(t) and the macroscopic current Jz(t)- We may summarize 
the relation as, 

J z (t) = J dt'a(t - t')E z (t') = --J dt'a(t - t')— jp, (24) 

where we have introduced the electric conductivity function cr(t). In the microscopic TDDFT 
calculation, the vector potential is an external variable. We may compute Jz(t) for an 
arbitrary Az(t) and thus determine the conductivity function. Since the equation is linear, 
it is easy to extract the frequency-dependent conductivity a{uj). The dielectric function e(u) 
is then given by the usual formula 

e{u) = 1 + (25) 
to 

Further details of calculating the dielectric function in this formalism are given in the Ap- 
pendix. 

Since our theory gives the same macroscopic current as in the linear response, the macro- 
scopic equations only require the dielectric function from the microscopic dynamics. Thus, 
the propagation of electromagnetic waves will be given by the usual relation between wave 
vector k and frequency u, 

u = * (26) 
te{u) 
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III. NUMERICAL 



A. Units 

For numerical quantities on the microscopic scale, we will use atomic units for length 
and field strengths. However, we will use micrometers for lengths on the macroscopic scale. 
On both scales, energies will be in eV units and time in femtosecond units. We continue to 
include the dimensionful quantities m and e in formulas even though their values are equal 
to one. For the laser intensity, we will use the conventional units of W/cm 2 . The conversion 
factor to atomic units is 1 a.u. = 3.509 x 10 16 W/cm 2 . 

B. Electronic scale 

As in our previous applications to crystalline materials (9), l^, we calculate the evolution 
of the electron wave function in a unit cell of the crystal. The orbital wave functions are 
represented on a 3D spatial grid which typically has a dimension of 16 3 . The Si lattice 
constant is 10.26 a.u. giving a mesh space of Ax = 0.64 a.u. A high order finite difference 
formula is used for the derivative calculations [35]. The number of fc-points in the reciprocal 
space cell is taken as 8 3 ; however due to symmetry there are only 80 distinct orbitals to be 
calculated. 

The number of ^-points adopted here are smaller than those employed in our previous 
work 0, [loj]. The present choice is decided from a computational feasibility. Our present 
scheme requires microscopic electron dynamics calculations in a number of macroscopic grid 
points simultaneously. Thus the present calculation consumes much more computational 
resources than our previous calculations of a single microscopic electron dynamics. The 
computational wall time with the above setting is approximately 15 hours employing 1,024 
processor cores in parallel with Intel Xeon X5570 (2.93GHz). This is close to the limit of 
our computational capability at present. The present calculation with 8 3 fc-points may not 
provide fully convergent results but we do not expect this truncation to affect the physical 
results by more than 10 percent. 

The electronic structure of each macroscopic position is initialized by the ground state 
Kohn-Sham orbitals. The wave functions are evolved by using the 4-th order expansion of 
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the TDDFT evolution operator |36|, |3 



e 



-iHM „ ^ (<At) B 



E K —^H\ (27) 



n=0 U - 



where H is the Kohn-Sham single-particle Hamiltonian appearing in Eq. f lT9|) . The orbitals 
at time t + At are computed by applying Eq. (1271) to the orbitals at time t. The Kohn-Sham 
operator if is a function of the fields 0, A, and the density. We use a fixed-time Hamiltonian 
H in which the scalar potential and the density are taken at time t and the field A is taken 
to be (A(t)+A(t + At))/2. This prescription does not require any predictor step (see below) 
and it gives good energy conservation over the course of the integration time. 
The algorithm (|2"T|) is stable provided At satisfies the condition 38] 

At < ^m(Ax) 2 « 0.2 a.u. (28) 

We use a somewhat smaller value in the calculations below, At = 0.08 au. In a typical run, 
the equations of motion are integrated for 16,000 time steps, amount to a total time 1280 au 
= 31 fs. This is sufficient to see the passage of a femtosecond laser pulse through a section of 
the solid. At longer times, other processes such as thermalizing collisions and ionic motion 
become important, and the TDDFT dynamics is no longer valid. 

As in previous work, we use the adiabatic approximation, taking the time-dependent 
functional in TDDFT the same as the ground state functional. We use the functional of 



Ref. 



in the local density approximation. The electron-ion interaction is treated with a 



norm-conserving pseudopotential 40|] with a gauge correction for the nonlocal part 



C. Macroscopic Scale 

On the macroscopic scale, Az(t) and J z (t) are considered as continuous functions, but 
they are discretized for the numerical calculation. In the results presented in the next section, 
we use a mesh size of 250 au ~ 13 nm. This permits us to propagate the pulse over a distance 
of several /im in the medium sampling the microscopic dynamics at several hundred points. 
We employ 256 grid points. The integrator for Eq. ( fl8~]) is straightforward, but the coupling 
between scales requires some care. The following update procedure is simple and conserves 
energy to adequate precision: 

f d 2 4vre 2 1 
A z (t + At) := 2A z (t) - A z {t - At) + c 2 At 2 -~Mt) + J z (t) , (29) 



dZ 2 
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where the space derivative is treated with a simple three-point formula. It also permits us 
to use Az(t + At) when updating the variables for the microscopic scale. 

D. Laser field 

We use the following functional form for the shape of the incident laser pulse, 

E z (t) = E sm( sm ), (Z < Z - ct < Z + cT). (30) 

cl c 

Here T = 18 fs controls the pulse width, ui = 1.55 eV/ft is the laser frequency, and E is the 
maximum electric field strength which is related to the laser intensity Iq by Iq = cEq/8tt. 

The corresponding gauge field is obtained by an analytic calculation of the following 
integral, 

A z (t) = -c f dt'E z {t'). (31) 

To start calculation, we need the initial vector potential at two times. One is given by 
Az{t = 0). Instead of using analytic form, we employ the following for the other, 

f)A 1 r) 2 A 

A z (At) = A z (0) + Ai^(0) + ^At^(0). (32) 

IV. RESULTS 

A. Pulse propagation 

We first note that the calculated dielectric constant at the laser frequency, e(u^) = 16.2, 
is in reasonable agreement with the observed value, e(ue) = 13.6. See the Appendix for 
details of the calculated dielectric function. The most significant shortcoming of the TDDFT 
dielectric function is this too-small band gap. Apart from that, we can confident that the 
present calculations will be reliable in the weak field limit. 

Snapshots of the time evolution for a typical run are shown in Fig. [TJ The initial laser pulse 
at t = has a peak intensity of Iq = 10 11 W/cm 2 and started at a position Z = —2.9 /im with 
respect to the Si surface at Z = 0. This is shown in the upper panel of the figure. The middle 
panel shows the field when the center of the pulse has just reached the surface, at t — 9.6 fs. 
One can see a transmitted wave of much smaller amplitude. In the lower panel, at t — 21.3 
fs, the wave has completely separated into the reflected and transmitted components. The 
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FIG. 1: Snapshots of the electromagnetic fields (vector potential divided by light speed, A/c, left 
panels) and of the electronic excitation energy per atom (right panels) at different times, shown as 
a function of macroscopic position. The vacuum is at Z < and the Si crystal is at Z > 0. Top 
panels: initial starting field, with pulse on left moving toward the Si surface. Middle panel: at 
the point the middle of the pulse reaches the surface. Lower panels: the reflected and transmitted 
pulses are well separated. The maximum intensity of the incident laser pulse is set 10 11 W/cm 2 . 

wave length of the transmitted component can be read off as A m = 3770 au, consistent 
with the low-field formula A m = A/ ' yfe ~ 3800 au. The center of the transmitted pulse is 
at Z = 0.71 fim. Taking the propagation time from the surface to be t2 — t\, the wave 
speed from the calculation is 0.20c. This is somewhat less than the phase velocity, which is 
c/y/e ~ 0.25c, but is completely consistent with the low-field group velocity computed as 



v> = v^m- l33) 

We also observe a chirp effect on the transmitted wave, stretched out at the front and 
condensed at the end of the transmitted pulse. 

We next examine the reflected and transmitted intensities. The maximum amplitude 
in Fig. [1] for the initial pulse is A /c = 0.0298, for the reflected pulse is A r /c = 0.0180, 
and for the transmitted pulse is A t /c = 0.0107. We obtain for the calculated reflectivity 
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r = (A r /Ao) 2 m 0.36. The reflectivity according to dielectric theory is given by 



at normal incidence. With our theoretical value for e(ui), we obtain R = 0.36, in good 
agreement with the real-time dynamics. The transmitted intensity is more complicated, 
since there are contributions from both the electronic part and the field part and the wave 
velocity is different. We can still ask how well the observed field amplitude agrees with 
dielectric theory. Expressing the transmittance T in terms of the field amplitudes, the 
formula is 



This gives T = 0.52 for the case shown in Fig. [U The dielectric transmittance can also be 
expressed purely in terms of e, giving T — 1 — R ~ 0.64. The difference between the two 
numbers, 0.64 — 0.52, is due to absorption. Thus the theory predicts that 12% of the energy 
is absorbed in the first 20 fs for an pulse of strength 10 11 W/cm 2 . In fact, as may be seen in 
the bottom right panel of Fig. [H we find a certain fraction of the excitation energy is left in 
the spatial region where the laser pulse already passed. Notice that this energy loss is not 
evident from the reflectance, which is still consistent with dielectric theory. 

In Fig. [21 we show energies per unit area integrated over the macroscopic coordinate. In 
the upper panel, the energy is decomposed into vacuum region (Z < 0, green dotted line) 
and Si crystal region (Z > 0, blue dashed line). The sum of the two contributions is shown 
by red solid line, showing that the total energy is well conserved during the whole period. 

In the lower panel, the energy per unit area in the Si crystal region is decomposed into 
contributions of electronic excitations and electromagnetic fields. Since the electromagnetic 
fields are separated into reflected and transmitted fields after 15 fs, the energy of Si crystal 
region does not change in that period. The energy of transmitted electromagnetic fields 
decreases gradually as it is transferred to electronic excitation. 

We next show reflected and transmitted electromagnetic fields at different intensity levels. 
In the left panels of Fig. [31 the vector potentials are shown at a time when the transmitted 
and reflected waves are well separated. In the right panels, the electronic excitation energies 
per atom is shown in the Si crystal region. At the lowest intensity, the propagation of electro- 
magnetic fields are well described by dielectric response. Essentially all of the energy remains 
associated with the propagating transmitted pulse. As the incident intensity increases, the 
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FIG. 2: Energies per unit area integrated over macroscopic coordinate Z are shown as a function 
of time. In the upper panel, the energies integrated over Z < (vacuum), Z > (Si), and the 
whole region (Total) are compared. In the lower panel, the energy integrated over Si crystal region 
is decomposed into the field part and the electronic excitation part. The incident laser pulse is the 
same as that of Fig. [TJ 
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FIG. 3: State of the system at t = 21 fs after the peak of the incident pulse reaches the surface for 
several different intensities of the incident laser pulse. Left: the field divided by light speed, A/c; 
right: excitation energy per atom in the Si crystal. 
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FIG. 4: The reflectivity of Si at normal incidence is shown as a function of peak laser intensity. 

transmitted wave becomes weaker than that expected from the linear response. We also find 
the central part of the transmitted pulse is suppressed strongly, producing a flat envelope 
of the pulse. In contrast, the envelope of the reflected wave does not change much in shape 
even at the highest intensity. We also find, at the intensity of 10 13 W/cm 2 , an emission 
of electromagnetic field is seen from the surface following the main pulse of reflected wave. 
From the right panels, above 10 12 W/cm 2 one sees that most of the energy is deposited in 
the medium with just a small fraction remaining in the transmitted electromagnetic pulse. 
The deposition rate falls off with depth as to be expected from the weakening of the pulse. 
At higher intensities the absorption rate greatly increases. At Iq = 10 13 W/cm 2 and higher 
the transmitted pulse is almost completely absorbed in the first tenths of a /xm. 

In Fig. HI we show the reflectivity as a function of incident laser intensity. Below 10 12 
W/cm 2 , the reflectivity is constant and in accord with dielectric theory (Eq. ( )34l) ). Above 
10 12 W/cm 2 , the reflectivity dips slightly, showing a minimum around 10 13 W/cm 2 . Above 
that intensity, the reflectivity start to increase gradually and finally reach 0.75 at the in- 
tensity of 5 x 10 14 W/cm 2 . This behavior of reflectivity qualitatively follows the observed 



evolution with intensity 20| , where it was interpreted in a dielectric model including effects 
of the excited electrons. We will later compare this model with our calculated reflectivity 
function. 
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FIG. 5: The vector potential divided by light speed (left panels) and electronic excitation energy 
(right panels) at the surface cell are shown as a function of time. 

B. Excitation in surface layer 

We next examine in more detail the first cell at the surface. The left-hand panel of 
Fig. [5] shows the vector potential as a function of time for several laser intensities. From a 
dielectric response, we expect the field inside the Si crystal is related to the incident field by 
At = (2/1 + y/e)Ai. This relation holds well below 10 12 W/cm 2 . At higher intensities, the 
field is less than this estimate gives. We also observe an oscillation of the vector potential 
after the incident pulse ends at 10 13 W/cm 2 , in accordance with what we found in Fig. [3j 
We will later consider this phenomenon with a model dielectric function. 

The electronic excitation energy in the first cell is shown in the right-hand panel of Fig. [51 
At the lowest intensities, the electronic energy is carried by the transmitted wave and leaves 
the cell after passage of the pulse. As the laser intensity increases, energy is transferred 
irreversibly to electronic excitation, and reaches a plateau after the laser pulse has passed 
(t > 15 fs). This is because the only mechanism to transfer energy between macroscopic 
grid points is through the macroscopic electromagnetic fields. 

Figure [6] shows some final-state properties of the surface as a function of intensity. The 
residual excitation energy is shown in the left-hand panel. At low intensities, the energy 
deposited is proportional to Jq. This is the expected dependence for two-photon absorption. 
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FIG. 6: Excitation of the first layer of Si after the laser pulse ends. Excitation energy per Si atom 
(left panel) and density of electron-hole pairs as number of pair per Silicon atom (right panel) are 
shown as a function of laser intensity. 

This is the most favorable absorption process in view of the photon energy: single-photon 
absorption is forbidden below the direct band gap, but the two-photon process is allowed. 
At Iq ~ 10 13 W/cm 2 the excitation energy is 0.6 eV per Si atom. This energy is in the 
form of electron-hole pairs. The minimum energy of a pair is at the direct band gap, 2.4 
eV. However, the excitation process forms a coherent pair with energies distributed across 
the valence and conduction bands. In the TDDFT dynamics, the coherence is lost after the 
pulse moves on, but the energy distribution remains the same. 

The number of particle-hole pairs n p h in the cell does not change after the electromagnetic 
field has passed. Then the number be calculated as the sum of overlaps of the time-dependent 
orbitals and the original Kohn-Sham orbitals, 



n 



* = £ i-Eife(o)fe(t)>f 



(36) 



where the sum over i,j is taken over occupied orbitals. The results are shown in the right 
panel of Fig. [6j As seen from the figure, the density increases quadratically with I Q up to 
to a point and then continues to increase more gradually. The ratio of energy density to 
particle-hole pair density, shown in Fig. [7J has a simple interpretation. At low intensities, 
up to about 10 12 W/cm 2 , it coincides accurately with two-photon energy 2hue = 3.1 eV. 
The energy per pair gradually increases at higher intensities. There one may expect two 
processes which increase the e nerg y per pair. One is higher-order multiphoton absorption, 



as has been often discussed 



41 



42]. The other is the secondary excitation of electrons which 
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FIG. 7: Electronic excitation energy per electron-hole pair as a function of laser intensity. 

have already been excited. 

With the information about the particle- hole density n p h, we may interpret the reflectivity 



curve (Fig. 4) with a mode 



example, in Ref. 



2fJ and 



for the dielectric function that includes plasma effects. For 



32|, the response of electrons excited in conduction band is 



described with the Drude model. We consider the following simplified form for the dielectric 
function, 

4:ire 2 n ph 1 



e(u,n ph ) = e(u,0) 



(37) 



Here e(u, 0) is the dielectric function in the ground state; m* and r are parameters of the 
Drude model. For our comparison we take e(u, 0) from the linear response (see Appendix) 
at uj = U£ = 1.55 eV. The reflectivity associated with the model dielectric function is 
determined from Eq. (134j) . Figure M shows the comparison for two assumptions about the 
effective mass and Drude damping time. For a given laser intensity, we use the electron-hole 
density n p h in our calculation shown in Fig. [6j The red open circle with solid line is the 
present calculation. The green filled circle with dotted line is the effective mass and damping 



0.18m and r = 1 fs. The blue filled circle with dashed 



time adopted in Ref. [20|, m* 
line is the parameters adopted in Ref. 32], m* = 0.35m and r = 0.5 fs. One can see that 
on a qualitative level, both the dip and the strong increase can be explained by plasma 
effects. One could try to fit the plasma parameters to reproduce the reflectivity curve, but 



19 
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FIG. 8: The reflectivity of Si at normal incidence is shown as a function of peak laser intensity. The 
red open circles with solid line repeat the calculated results from Fig. 0J The green filled circles 
with dotted line and blue filled circles with dashed line use Eqs. (|34|) and (137p with m* = 0.18m, 
r = l.Ofs, and m* = 0.35m, r = 0.5fs, respectively. 

it is probably not realistic to assume that a fixed dielectric function is responsible for the 
electromagnetic interactions. However, it should be mentioned that the reflectivity as well 
as the absolute value of the dielectric function are minimized when the screened plasma 
frequency, 

2 _ 4Tie 2 n ph 
~ e(u;„0)m*' (38) 

coincides with the frequency of the incident laser pulse, = u p . This relation is fulfilled at 
the laser intensity around 10 13 W/cm 2 , consistent with the behavior of reflectivity. 

In Figs. [3] and [3, we observed an emission of electromagnetic field following the main 
pulse of the reflected wave at the laser intensity of 10 13 W/cm 2 . This phenomenon may 
also be understood with the model dielectric function. At this intensity, a small magnitude 
of the dielectric function at the surface allows a penetration of transmitted wave inside 
the medium. However, the dielectric function changes rapidly inside the medium due to 
the increase of electron-hole pair density. It may cause a reflection from a deeper layer, 
producing the electromagnetic field following the main pulse. 



20 



10" 



I 10° 
> 



^o'■ 



1— 

CD 
c 

CD 
"D 

B 10" 4 
o 

X 

LU 



10'* 




Multi-Scale — e- 
Longitudinal -----©-- 
Transverse o~ 



10 a 



10 10 10 11 



10 12 10 13 



1Q 14 1Q 15 1Q 16 



Intensity [W/cm2] 



10" 



I 10° 
> 

CO 



10"' 



1— 
CP 

c 

CD 
"D 

B 10" 4 
o 

X 

LU 



10'* 



10 



10 



Multi-Scale — e- 

Longitudinal &-■ 

Transverse 
E=Cxl 2 



10' 



10 



12 



10 



13 



10 



14 



Intensity [W/cm2] 



10 



15 



FIG. 9: Deposited energy in the Si medium. Red solid line: the energy deposited in the first- 
layer in the multi-scale calculation. Green, dashed line: microscopic calculation in the longitudinal 
geometry. Blue dotted line: microscopic calculation in the transverse geometry. In the left panel, 
the horizontal axis is the intensity of the incident laser pulse for multi-scale calculation, and is the 
intensity of the applied laser pulses in the microscopic calculations of longitudinal and transverse 
geometries. In the right panel, the laser intensity is normalized to the transverse case. See the text 
for more detail. 



C. Multi-scale vs single-cell approximations 

In Ref. {9 1, we calculated microscopic electron dynamics for an external electric field 
normal to the crystal surface and neglecting magnetic fields. In this longitudinal geometry, 
the crystal response is uniform and the TDDFT is computationally much less expensive. This 
was applied to the dielectric breakdown for diamond crystal, and the calculated threshold 
for breakdown was at least an order of magnitude higher than the measured threshold. 

The aim of the present subsection is twofold. First we show that the present multi-scale 
calculation gives much lower breakdown threshold than that of our previous calculation 
in the longitudinal geometry, thus resolving the discrepancy of our previous calculations 
with measurements. Second, we clarify mutual relationship between the present multi-scale 
calculation and the single-cell treatment in either the longitudinal or transverse geometry. 
Since the single-cell calculations are much easier computationally, it would be useful to know 
what physical information may be extracted reliably from them. 

We first explain in more detail the longitudinal and transverse geometries in the single- 
cell calculations. In the transverse geometry, we simply put the vector potential of applied 
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laser pulse, A(t), in the Kohn-Sham Hamiltonian and calculate the electron dynamics. In 
the longitudinal geometry, we take that field as external and add to it the field from the 
induced current in the medium. The vector potential in the Kohn-Sham Hamiltonian is the 
sum of the external and the induced fields, A(t) = A ext (t) + A in( i(t). 

The final-state electronic excitation energies for the three calculations are shown in Fig. 
[9j In the left panel, the red circles and solid line shows the deposited energy at the surface in 
the multi-scale calculation as a function of the incident laser intensity. The green circles and 
dashed line is the microscopic calculation in the longitudinal geometry, as adopted in Ref. 
9|. The blue circles and dotted line is the microscopic calculation in the transverse geometry. 
We may identify the dielectric breakdown at the laser intensity where the electron excitation 
energy per atom is about 1 eV. One sees that the threshold for dielectric breakdown is very 
different for the three calculations. The threshold is lower by an order of magnitude for the 
multi-scale and transverse cases than the longitudinal case. 

The difference may be understood using a dielectric picture to relate the internal and 
external fields. In the transverse case, the applied electric field directly acts upon electrons 
in the medium. In the case of the multi-scale calculation, the electric field in the medium 
and the incident field are related by 

2 

^medium Z j ~~F-£in (^9) 
-L + V 6 

Putting the value of dielectric constant e = 16, the laser intensity is different between the 
transverse and multi-scale calculations by a factor of (2/5) 2 = 0.16. In the longitudinal case, 
in addition to the above factor connecting medium and incident fields, we need to add the 
following factor connecting the external and the medium fields, 

£-ext ^medium j (40) 

The factor to correct the laser intensity is 16 2 (2/5) 2 = 41 for the longitudinal geometry. 
Taking these factors as corrections to the laser intensity, we replot the electronic excita- 
tion energy as a function of laser intensity in the medium in the right panel of Fig. 
We see that these factors explain accurately the order-or-magnitude difference in the di- 
electric breakdown threshold. The electronic excitation energy coincides accurately below 
10 12 W/cm 2 for three calculations, where excitations are mostly by two-photon absorption. 
There are some deviations around 10 13 W/cm 2 and above, where the resonant excitation is 
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expected. The longitudinal calculation shows an abrupt rise of the excitation energy which 
we interpreted as a resonant energy transfer from the laser pulse to the electrons [9|. The 
other two calculations do not produce an abrupt rise but rather show a smooth saturation 
of the energy transfer. 



V. SUMMARY AND OUTLOOK 



We have developed a first-principles framework to calculation the propagation of electro- 
magnetic field in crystalline solids. The macroscopic electromagnetic field is described by 
Maxwell equations while the microscopic electron dynamics is described by TDDFT. With 
use of massively parallel computers, we showed that it is feasible to treat one of the simplest 
systems of physical interest, the propagation of a laser pulse into bulk Si at the normal 
incidence. 

At low field intensity, the calculated field propagation and electronic excitations exhibit 
features expected from ordinary electromagnetic theory with the dielectric function given by 
linear response theory. The electronic excitations are dominated by two-photon absorption 
at low intensities since the laser frequency is below the direct bandgap. 

As the laser intensity increased, the density of excited electron-hole pairs become high 
enough to affect the response. This is conveniently modeled as an electron-hole plasma. At 
around 10 13 W/cm 2 , the plasma frequency of excited electrons reaches the visible frequency, 
showing a nonlinear interaction with the incident laser pulse. Above this intensity, the 
responses are dominated by nonlinear electron dynamics. 

We have also found that the surface absorption obtained in the multi-scale theory can be 
described by a single-cell approximation using dielectric formulas to relate the internal and 
external fields, provided the fields do not much exceed 10 13 W/cm 2 . 

Finally, we mention some directions that might be interesting to take up in later work. 
Analytic approximations have been proposed to express the excitation energy as a function 
of the Keldysh parameter [41]]. We have not examined the validity or accuracy of such 
approximations, but it would be useful to have this information. 

Computations in the present framework could be extended to deal with laser pulses at 
oblique angles of incidence. In that case, the field is not translationally invariant in the x 
direction, but the medium itself is. Consequently relatively few cells would be needed to 
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describe the x dependence. It would also be interesting to extend the present calculations to 
pump-probe laser pulse protocols. In principle it is straightforward to calculate the response 
to a double pulse separated in time. As a practical matter, pump-probe responses could 
most easily be studied in a single-cell approximation. Also, one could examine the linear 
response of the excited system using fields of the pump-probe form. This is important to 
verify the validity of the arguments made in Sect. IIVBI 
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Appendix 

Here we show how the dielectric function may be calculated using the formalism of Sect. 
II. E. For the perturbation, we take A to be of the form 

A{t) = A 9(t) . (41) 

The microscopic equation of motion Eq. (Q is integrated from t — to t — T m to obtain the 
J(t) over the time interval. As is evident from Eq. f)24p . the calculated J(t) is proportional 
to the conductivity as a function of time, 

a(t) = -j-J(t). (42) 

This is Fourier transformed as 

fTm 

a{u) = / dte lujt f(t)a(t) (43) 



o 



where f(t) is a filter to suppress spurious oscillations that would arise 



10m a sharp cutoff 



of the integration at T m . We employ a third order polynomial for it 3^|. The dielectric 
function may be obtained from the conductivity by Eq. (125]) . 
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FIG. 10: Conductivity as a function of time. Calculations for two choices of fc-points are compared. 
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FIG. 11: Conductivity and dielectric function as a function of frequency in which 32 3 fc-points are 
used. The measured value is also shown for dielectric function. 

We carried out this computation taking A Q = 0.0005 a.u. and T m = 31 fs (16,000 time 
steps with At = 0.08 a.u.). In Fig. [101 we show a conductivity as a function of time, cr(t), 
for two choices of fc-points, 8 3 and 32 3 . In our multi-scale calculation, we adopt 8 3 fc-points. 
Two calculations coincide each other up to 2 fs. There remain oscillations for a long period 
in the calculation of 8 3 fc-points, which are washed out if one employs a finer fc-points grid. 
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The conductivity a(t) is Fourier transformed to obtain the conductivity and dielectric 
function as a function of frequency. They are shown in Fig. [Ill i n which 32 3 fc-points are 
used. At a frequency region close to zero, the conductivity ciui) should behave 

da u. (44) 



a[U) = 1 fa 



w=0 



In actual calculation, it is not exact due to the presence of spurious mode which originates 
from a violation of translational invariance in the real-space grid calculation. Since a small 
deviation from the above analytic behavior at around u = harms the low frequency 
behavior of the dielectric function, we replace the real part of the dielectric function by a 
second order polynomial of the frequency below 1 eV. The calculated dielectric function, 
e(u), shown in the right panels of Fig. [TTJ is very close to the one calculated in Ref. [uj] 



using the formalism of Ref. 



34]. 



The calculated real part of the static dielectric function is e(0) = 12.6, close to the 
experimental value of 11.6. However, as is well known in the density functional theory, the 
direct band gap in the local density approximation is smaller than the experimental one (2.4 
eV theory vs. 3.1 eV experiment). 
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